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SPACE-TIME  CHARACTERISTICS  OF  THE  TURBULENT 
FLOW  ENERGY  IN  A  VISCOELASTIC  MEDIUM 
USING  NONINVASIVE  PASSIVE  IMAGING 

CHAPTER  1.  INTRODUCTION 

Attempts  to  noninvasively  diagnose  coronary  artery  disease  (Akay  et  al.,  1992)  and 
intracranial  vascular  lesions  and  aneurysms  (Abiko  et  al.,  1994)  by  passive  listening  (auscultation) 
for  the  effects  of  blood  flow  turbulence  have  been  reported  in  the  recent  literature.  The  potential 
for  these  techniques  to  result  in  low  cost,  mass  screening  examinations  is  considerable.  However, 
the  feasibility  of  such  noninvasive  detection  of  human  arterial  disease  for  other  than  a  near-surface 
artery  depends  on  the  ability  to  discriminate  signals  with  a  very  low  energy  level  from  high  level 
background  masking  interference  and  physiologically  produced  noise.  It  is  proposed  here  to 
enhance  the  desired  signal-to-noise  ratio  (SNR)  through  the  use  of  an  array  of  noninvasive 
vibroacoustic  sensors  placed  on  the  body  surface.  It  is  suggested  that  the  wave  lengths  for 
vibrations,  which  occur  in  the  frequency  decade  from  100  to  1000  Hz  (Verberg,  1983),  are 
appropriate  for  spatially  discriminating  signals  from  noise  on  the  basis  of  joint  space-time  signal 
processing.  The  space-time  signals  of  interest  result  from  (blood)  turbulence-induced  forces 
normal  and  parallel  to  a  diseased  artery  internal  wall,  which,  in  turn,  couple  as  wave  energy  to  the 
body  surface  by  means  of  a  time-varying  spatial  vibration  field.  This  space-time  field  describes  the 
energy  propagation  in  the  viscoelastic  body  medium.  It  is  proposed  that  the  three-dimensional 
vibration  field  can  be  sensed  by  an  array  of  transducers  placed  on  the  skin.  The  spatial 
discrimination  properties  of  the  array  would  be  used  as  a  generalized  focusing  "antenna"  to 
enhance  the  SNR  relative  to  that  which  could  be  achieved  with  a  single  sensor  (Owsley  et  al., 
1985;  Owsley,  1994). 

A  blood  flow  turbulence-induced  murmur  can  be  discriminated  from  masking  background 
sound  in  the  temporal,  spectral,  and  spatial  domains.  For  example,  temporal  discrimination  against 
background  noise  in  auscultation  for  coronary  artery  disease  (CAD)  is  based  on  the  simple  fact  that 
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sounds  characterizing  occlusive  coronary  stenosis  occur  during  diastole  (Akay  et  al.,  1990). 
Energy  from  stenosed  coronary  arteries  must  therefore  be  extracted  from  a  background  noise 
environment  that  results  from  the  rapid  filling  of  the  ventricle  chambers,  including  the  first  and 
second  heart  sounds.  Moreover,  it  has  been  suggested  that  resonance  effects  in  arteries  excited  by 
blood  flow  produce  characterizable,  albeit  low  level,  sonic  signatures  (Wang,  1990).  Double 
resonant  frequencies  below  400  Hz  are  predicted,  whereas  experimentally  CAD-related  resonant 
energy  as  high  as  1 100  Hz  due  to  turbulent  blood  flow  caused  by  lesions  on  the  inner  wall  of  an 
occluded  artery  is  reported  (Akay  et  al.,  1992).  In  addition,  quantitative  blood  flow  rate  may  be 
measurable  in  the  absence  of  highly  turbulent  flow-induced  resonance  by  estimating  the  level  and 
cutoff  frequencies  of  the  broadband  energy  produced  by  low  level  turbulence,  which  is,  in  turn, 
related  to  the  diameter  of  the  stenotic  carotid  segment  (Lees  and  Kistler,  1978).  In  this  study,  it  is 
suggested  that  a  key  to  reliable  detection,  identification,  and  location  of  arterial  lesions  and 
aneurysms  may  very  well  lie  in  the  spatial  domain.  It  is  thus  proposed  that  the  detection  of  the  type 
and  location  of  both  normal  and  pathologic  blood  flow  vibration-induced  acoustic  and  nonacoustic 
energy  can  potentially  be  accomplished  through  simultaneous  filtering  and  combining  of  the 
outputs  of  a  spatially  diverse  array  of  vibro- acoustic  transducers. 

For  CAD,  the  signal  levels  produced  by  the  effects  of  relatively  small  coronary  artery  blood 
flow  rates  are  substantially  below  the  levels  produced  by  nearby  flow  activity,  such  as  ventricle 
refilling  during  diastole.  Moreover,  the  ambient  environmental  noise  produced  by  noncardiac 
functions  is  a  known  inhibitor  to  low  level  sound  detection  and  diagnosis.  It  is  well  known  that 
the  audible  sonic  field  produced  by  the  cardiac  cycle  exhibits  a  high  degree  of  spatial  variability 
(Epstein,  1991).  The  traditional  auscultation  in  phonocardiography  involves  the  placement  of  a 
single  vibro-acoustic  transducer  anywhere  in  the  area  of  energy  propagation  associated  with  high 
quality  transmission  from  a  particular  region  of  the  heart.  In  contrast,  multiple  auscultation  point 
(MAP)  space-time  phonocardiography  proposes  to  receive  energy  simultaneously  at  an  array  of 
transducers  and  through  the  application  of  modem  space-time  signal  processing  technology  to 
provide: 
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(a)  a  maximum  SNR  estimate  of  the  waveform  produced  by  the  turbulent  blood  flow 
activity  in  a  specific  region  of  the  body  by  coherently  combining  the  spatially  diverse  samples  of 
the  vibration  field,  thereby  enhancing  the  desired  signal  and  discriminating  against  interfering 
components  to  the  measured  space-time  field  and 

(b)  a  spatially  and  spectrally  resolved  image  of  the  turbulent-flow-induced  spatially 
extended  vibratory  energy  source  for  either  a  narrowband  (resonant)  phenomenon  or  a  spectrally 
wideband  phenomenon  (e.g.,  low  intensity  level  turbulent  flow  activity  and  artifacts  such  as 
stenotic  valve  clicks  and  prosthetic  valve/vessel  failure  points. 

Enhanced  SNR  waveform  estimation  is  critical  to  low  error  rate  detection  algorithms  for 
low  level  sonic  activity.  High  resolution  imaging  of  energy  fields  due  to  turbulent  fluid  jets  would 
be  invaluable  for  the  prognosis  and  treatment  of  lesions  and  aneurysms.  The  combination  of  these 
two  capabilities  has  the  potential  for  a  relatively  low  cost  alternative  to  invasive  procedures  and  in 
particular  to  catheterization. 

This  report  first  defines  the  experimental  apparatus  and  array  transducer  used  to  measure 
the  viscoelastic  properties  that  support  wave  propagation  in  the  body  resulting  from  turbulent-flow- 
induced  vibrations.  Next,  a  proof-of-concept  (POC)  experiment  that  establishes  the  feasibility  of 
the  proposed  MAP  concept  using  frequency- wavenumber  analysis  is  described,  and  an  attempt  at 
the  interpretation  of  the  experimental  results  in  terms  of  propagating  shear  and  surface  waves  is 
examined.  The  potential  for  three-dimensional  high  resolution  imaging  is  validated  through  a 
numerical  simulation  of  the  simulated  space-time  vibrational  energy  field.  Finally,  work  is 
proposed  that  constitutes  the  basis  for  the  development  of  a  commercially  viable  clinical  application 
of  the  MAP  concept  to  very  low  cost,  noninvasive  mass  diagnostic  screening. 
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CHAPTER  2.  EXPERIMENTAL  APPARATUS 

2.1.  LABORATORY  CONFIGURATION 

The  experimental  apparatus  consists  of  a  10.9  by  6.1  by  4.0  inch  rectangular  Hexcel  195- 
RE  urethane  block  with  embedded  surgical  tubing  at  a  depth  of  0.69  inch  (figure  2. 1).  Inside  the 
surgical  tubing  is  a  cylindrical  nylon  blockage  that  is  1.0  inch  long  with  two  U-shaped  cuts,  as 
shown  in  figure  2.2.  The  outer  diameter  of  the  blockage  is  0.125  inch  and  the  inner  diameter  is 
0.043  inch.  The  blockage  is  designed  to  convert  an  incoming  laminar  flow  field  into  a  turbulent 
flow  field.  An  elevated  50-gallon  upper  reservoir  feed  tank  is  located  next  to  the  urethane  block 
and  connected  to  the  surgical  tubing.  The  water  drains  from  the  upper  reservoir  through  the 
urethane  block  and  into  a  calibrated  flask  located  below  the  block  (figure  2.3).  The  fluid  velocity 
can  be  changed  by  varying  the  head  pressure  of  the  50-gallon  tank.  One  sensor  was  positioned  at  a 
fixed  location  directly  above  the  blockage  and  the  other  was  placed  at  various  locations 
downstream  of  the  blockage  during  the  course  of  the  experiment. 

2.2.  SENSOR 

The  sensor  used  in  this  experiment  is  a  polyvinylidene  fluoride  (PVDF)  motion  sensor 
(Kassal  et  al.,  1994).  The  PVDF  piezoelectric  polymer  produces  charges  of  equal  magnitude  and 
opposite  polarity  on  opposite  surfaces  when  a  mechanical  strain  is  imposed  on  the  material.  The 
size  of  the  surface  charge  is  proportional  to  the  size  of  the  strain.  The  sensor  is  designed  for  use  in 
cardiac  auscultation,  although  it  measures  motion  rather  than  acoustic  pressure.  It  is  1  inch  by 
1.25  inches  and  the  active  (measurement)  size  of  the  sensor  is  1  inch  by  1  inch.  It  consists  of  two 
rectangular  thin-film  PVDF  elements,  a  flexible  substrate,  a  miniature  three-conductor  electrical 
connector,  a  moisture  barrier  that  also  acts  as  a  protective  coating,  and  a  connector  bib  to  keep  the 
cable  connector  receptacles  clean.  Hydrogel  adhesive  is  applied  between  the  sensor  and  the 
measurement  location  to  ensure  good  mechanical  contact.  The  signal  from  the  sensor  is  passed  to  a 
differential  amplifier  to  increase  the  signal  and  then  into  an  HP3562A  signal  analyzer. 


4 


TR  10,833 


Figure  2. 1.  Urethane  Block  With  Sensors 
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CHAPTER  3.  EXPERIMENTAL  RESULTS  AS  AN 
INITIAL  PROOF  OF  CONCEPT  (POC) 

3.1.  INTRODUCTION 

An  initial  proof-of-concept  experiment  was  performed  on  the  laboratory  apparatus 
described  in  chapter  2.  The  results  were  obtained  by  first  measuring  frequency  domain  transfer 
functions  between  a  stationary  sensor  and  a  relocatable  sensor  positioned  at  15  locations  and  then 
by  applying  a  Fourier  transform  to  the  spatial  domain  data  for  conversion  into  the  wavenumber 
domain.  This  process  is  equivalent  to  that  of  a  farfield  beamformer  whose  focal  point  is  at  infinity. 
These  frequency-wavenumber  (kco)  data  surfaces  are  displayed  primarily  as  wavenumber  cuts  at 
constant  frequency.  Wave  motion  is  identified  on  these  wavenumber  cuts,  and  the  corresponding 
apparent  wave  speed  is  calculated.  These  wave  speeds  are  then  compared  to  theoretical  values 
determined  from  material  property  measurements.  The  energy  levels  present  in  the  sidelobe 
structures  of  the  beam  are  also  analyzed. 

3.2.  MATERIAL  MEASUREMENTS 

The  Hexcel  195-RE  urethane  and  the  surgical  tubing  were  tested  using  the  Metravib 
Viscoanalyseur  in  order  to  determine  their  dynamic  material  properties.  The  Viscoanalyseur  is  a 
material  testing  machine  designed  to  measure  the  complex  compressional  and  shear  moduli  as  a 
function  of  frequency  and  temperature  (Metravib,  1989).  Table  3.1  is  a  listing  of  the 
compressional  modulus  and  table  3.2  is  a  listing  of  the  shear  modulus  of  Hexcel  195-RE  urethane. 
Table  3.3  lists  the  compressional  modulus  and  table  3.4  lists  the  shear  modulus  of  the  surgical 
tubing.  All  four  tables  show  measurements  made  at  20  degrees  C.  For  completeness,  we  note  that 
using  the  compressional  and  shear  moduli  to  calculate  a  value  for  Poisson's  ratio  yields  a  number 
ranging  from  2  to  4,  a  physically  unrealizable  value.  It  is  believed  that  this  error  resides  in  the 
measured  values  for  the  compressional  modulus.  The  Materials  Laboratory  at  NUWC  Detachment 
New  London  is  addressing  the  issue  of  incompatible  Poisson's  ratios.  The  density  of  Hexcel  195- 
RE  urethane  is  930  kg/m^. 
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Table  3.1.  Compressional  Modulus  of  Hexcel  195-RE  Urethane 
Measured  With  the  Viscoanalyseur 


Frequency  (Hz) 

5.0 

54.5 

104.0 

153.6 

202.9 

252.5 

301.9 
351.1 

400.6 
451.3 
500.0 


Real(E)  (N/m^)  Iraaginary(£)  (N/m^) 


4.90  X  10^ 
6.29  X  10^ 
8.26  X  10^ 
9.03  X  105 
9.68  X  105 

10.10  X  10^ 

10.50  X  10^ 
10.70  X  105 
10.80  X  10^ 
10.50  X  105 
9.72  X  10^ 


1.25  X  10^ 

3.80  X  10^ 

3.47  X  10^ 
3.95  X  10^ 
4.28  X  105 
4.59  X  10^ 
4.89  X  105 
5.17  X  lO^ 
5.39  X  10^ 
5.67  X  10^ 
6.30  X  105 


Table  3.2.  Shear  Modulus  of  Hexcel  195-RE  Urethane 
Measured  With  the  Viscoanalyseur 


Frequency  (Hz) 

Real(G)  (N/m^) 

Imaginary(G)  (N/m^) 

5.0 

1.04  X  10^ 

0.149  X  10^ 

54.5 

1.32  X  10^ 

0.313  X  10^ 

104.0 

1.43  X  10^ 

0.386  X  10^ 

153.6 

1.52  X  10^ 

0.438  X  10^ 

202.9 

1.59  X  105 

0.466  X  10^ 

252.5 

1.65  X  10^ 

0.505  X  10^ 

301.9 

1.69  X  10^ 

0.537  X  10^ 

351.1 

1.73  X  10^ 

0.557  X  105 

400.6 

1.67  X  105 

0.555  X  10^ 

451.3 

1.63  X  10^ 

0.618  X  10^ 

500.0 

1.45  X  10^ 

0.698  X  10^ 
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Table  3.3.  Compressional  Modulus  of  Surgical  Tubing 
Measured  With  the  Viscoanalyseur 


Frequency  (Hz) 

Real(£)  (N/m^) 

Iraaginary(£)  (N/m^) 

5.0 

69.4  X  10^ 

4.59  X  10^ 

54.5 

76.1  X  10^ 

6.03  X  10^ 

104.0 

77.3  X  10^ 

7.13  X  10^ 

153.6 

77.2  X  10^ 

9.21  X  10^ 

202.9 

80.2  X  10^ 

8.03  X  10^ 

252.5 

81.3  X  10^ 

8.29  X  105 

301.9 

82.4  X  10^ 

8.43  X  10^ 

351.1 

83.4  X  10^ 

8.72  X  105 

400.6 

84.3  X  10^ 

8.97  X  10^ 

451.3 

85.1  X  10^ 

9.26  X  10^ 

500.0 

85.8  X  10^ 

9.68  X  10^ 

Table  3.4.  Shear  Modulus  of  Surgical  Tubing 
Measured  With  the  Viscoanalyseur 


Frequency  (Hz) 

Real(G)  (N/m^) 

Imagmary(G)  (N/m^) 

5.0 

7.45  X  10^ 

0.576  X  105 

54.5 

7.93  X  10^ 

0.843  X  10^ 

104.0 

7.94  X  10^ 

1.010  X  10^ 

153.6 

8.02  X  10^ 

1.150  X  10^ 

202.9 

8.19  X  10^ 

1.160  X  10^ 

252.5 

8.29  X  10^ 

1.290  X  10^ 

301.9 

8.39  X  10^ 

1.360  X  10^ 

351.1 

8.43  X  10^ 

1.380  X  10^ 

400.6 

8.48  X  10^ 

1.500  X  10^ 

451.3 

8.45  X  10^ 

1.570  X  105 

500.0 

8.28  X  10^ 

1.690  X  105 
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Based  on  these  shear  modulus  measurements  note  that  the  shear  wave  speed  in  the  urethane 
is  10  to  14  m/s  in  the  frequency  region  of  0  to  500  Hz.  The  shear  wave  speed  of  typical  tissue  is 
13  m/s  at  1000  Hz  (Verburg,  1983)  although  the  dispersive  properties  of  the  shear  wave  in  tissue 
are  more  severe  (at  lower  frequencies)  than  the  urethane  used  in  this  test.  Because  the 
compressional  wave  speed  in  the  urethane  is  difficult  to  calculate  due  to  the  uncertainty  of  the 
measured  compressional  modulus,  the  measurements  of  the  POC  experiment  are  shown  (below)  to 
be  shear  wave  dominant  and  compressional  energy  to  be  negligible. 

3.3.  FREQUENCY-SPATIAL  DOMAIN  MEASUREMENTS 

The  frequency-spatial  domain  measurements  were  made  using  an  HP3562A  dynamic  signal 
analyzer.  The  analyzer  converts  the  raw  data  from  the  time  domain  to  the  frequency  domain.  The 
measurements  were  made  using  two  different  measurement  types:  power  spectral  density  and 
transfer  function.  Power  spectral  density  is  a  measure  of  the  energy  of  a  single  sensor,  and  the 
transfer  function  is  a  measure  of  the  relative  amplitude  and  phase  in  two  sensors.  Both 
measurements  are  made  with  respect  to  frequency.  We  define 


(3.1a) 

Gyy{f)=Fyif)Fyif)  , 

(3.1b) 

(3.1c) 

where  G^xif)  and  Gyy(f)  are  the  autopower  spectrum  densities  (APSD)  of  channel  x  and  y, 
respectively;  G^{f)  is  the  crosspower  spectrum  density  (CPSD);  /^(/)  and  /^(/)  are  the  linear 

^  9(C 

complex  spectra  of  the  waveforms  fy^{t)  and  fy{t),  respectively;  if)  and  Fy  (/)  are  the 

complex  conjugate  of  the  linear  spectrum; /is  the  frequency;  and  the  bar  above  the  right-hand  side 
denotes  a  time  average.  The  transfer  function  between  channels  x  and  y  is  defined  as 
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Hxy(f)  = 


Gxy(f) 

Gxxif) 


FxWFyif) 

Fxif)Fx{f)  ' 


(3.2) 


Additionally,  the  magnitude-squared  (MSC)  coherence  function  between  channels  x  and  y  is 
defined  as 


2  _  Gxyif)Gxyif) 

Gxx(/)%(/)  ■ 

Figure  3.1  is  a  plot  of  the  APSD  of  a  sensor  pad  directly  above  the  blockage.  The  solid 
hne  shows  data  taken  with  the  flow  and  the  dashed  line  shows  data  without  flow.  Note  that  at 
some  frequencies  there  is  almost  a  20-dB  SNR.  Figure  3.2  plots  the  transfer  function  between 
sensors  whose  phase  centers  are  separated  by  1.09  inches.  Because  the  distance  from  the  phase 
center  to  the  edge  of  each  sensor  is  approximately  0.5  inch,  figure  3.2  depicts  the  sensors' 
placement  as  being  almost  contiguous.  The  first  sensor  in  figure  3.2  was  located  directly  above  the 
blockage  and  the  second  sensor  was  located  downstream.  Although  the  noise  flow  defined  by  the 
no-flow  condition  is  not  shown  in  this  or  the  subsequent  transfer  functions,  it  is  at  an 
approximately  -15-dB  level  in  the  magnitude  plot  of  the  transfer  functions.  Figure  3.3  is  a  plot  of 
the  transfer  function  between  sensors  whose  phase  centers  are  separated  by  2.75  inches.  Figure 
3.4  is  a  plot  of  the  transfer  function  between  sensors  whose  phase  centers  are  separated  by  4.50 
inches.  In  figures  3.3  and  3.4,  the  second  sensor  was  moved  downstream  and  the  first  remained 
over  the  blockage.  The  phase  angle  plots  in  figures  3.2,  3.3,  and  3.4  show  well-defined 
propagating  energy  content.  Figures  3.5,  3.6,  and  3.7  are  plots  of  the  MSC  functions  with  sensor 
phase  centers  separated  by  1.09,  2.75,  and  4.50  inches,  respectively.  Although  the  coherence 
decreases  as  the  measurement  distance  from  the  blockage  increases,  the  signal  is  still  highly 
coherent  at  36  tube  diameters  downstream.  All  the  data  shown  here  were  taken  using  a  frequency 
range  of  0  to  1000  Hz.  The  500-  to  1000-  Hz  portion  of  the  data  is  not  shown  because  the  sensor 
is  spatially  averaging  the  high  frequency  content  and  coherent  energy  structures  are  removed.  A 
smaller  sensor  will  not  produce  this  effect.  For  all  of  the  measurements,  the  flow  rate  was 
0.00348  liters/second,  which  has  a  corresponding  Reynolds  number  of  1400.  Although  this  value 
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is  below  the  Reynolds  number  for  turbulent  flow  in  a  pipe,  the  flow  field  around  the  blockage  will 
be  turbulent  due  to  the  change  in  area  caused  by  the  blockage. 
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Phase  Centers  Separated  by  4.50  Inches 


3.4.  TRANSFORMATION  FROM  SPATIAL  TO  WAVENUMBER  DOMAIN 

The  wavenumber  domain  allows  the  display  of  data  that  are  dominated  by  several  waves. 
As  described  above,  the  data,  collected  at  each  sensor  channel  in  the  form  of  a  frequency  domain 
transfer  function  of  the  downstream  sensor  divided  by  the  stationary  sensor  above  the  blockage, 
are  transformed  from  the  spatial  domain  into  the  wavenumber  domain  by  means  of  a  discrete 
spatial  Fourier  transform  at  each  frequency.  This  is  the  same  process  that  a  conventional 
beamformer  (CBF)  uses  to  detect  propagating  farfield-source  planar  wave  energy. 

The  discrete  Fourier  transform  originates  with  the  Fourier  integral,  which  converts  a 
continuous  function  from  the  spatial  (x)  domain  (meters)  to  the  wavenumber  (k)  domain 
(radians/meter)  (Corcos,  1963;  Bradshaw,  1971;  Smol'yakov  and  Tkachenko,  1983): 

H(k)  =  —  J h(x)e~^^dx  ,  (3.4) 

2  TV 
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where  h{x)  is  the  function  in  the  spatial  domain,  H(k)  is  the  function  in  the  wavenumber  domain,  x 
is  distance,  k  is  wavenumber,  and  j  is  the  square  root  of  -1.  The  left-hand  side  of  equation  (3.4) 
has  units  that  are  equal  to  the  units  of  h(x)  multiplied  by  meters/radian.  The  limits  of  equation 
(3.4)  can  be  changed  to  reflect  a  finite  length  interval,  which  is  written  as 


,  L 

Hik)  =  —\hix)e~j^dx  , 
271  i 


(3.5) 


where  L  is  the  length  of  the  interval. 

In  cases  involving  multiple  discrete  sensors,  it  is  not  possible  to  sample  h{x)  continuously 
on  the  aperture  interval.  It  is  possible,  however,  to  transform  the  measurements  from  the  spatial 
domain  to  the  wavenumber  domain  by  means  of  a  discrete  Fourier  transform  having  the  form 
(Bendat  and  Piersol,  1980;  Bendat  and  Piersol,  1986) 


N-l 

Hi(k)  =  H(iAk)=-^J^hiXn)e 


(3.6) 


where  Hiik)  is  the  discrete  wavenumber  domain  transform,  h(x^)  is  the  discrete  spatial  domain  data 
at  the  nth  sensor,  (=  UN)  is  the  separation  distance  between  discrete  sensors  (m),  N  is  the  total 
number  of  sensors,  i  is  a  counter  that  runs  from  1  to  N,  and  Aj^  is  the  discrete  wavenumber 
spacing  of  the  transformed  data  (radians/meter).  The  wavenumber  spacing  is  related  to  the  aperture 
length,  the  spacing  of  the  sensors,  and  the  number  of  sensors  by 


1%  _  2% 

^~~L~  NA^ 


(3.7) 


Equation  (3.6)  provides  a  means  to  transform  data  sampled  with  multiple  sensors  into  the 
wavenumber  domain. 

A  process  known  as  zero  padding  can  be  used  to  interpolate  the  data  in  wavenumber.  This 
process  involves  using  an  additional  M  spatial  data  points  in  equation  (3.6),  all  of  which  are  set 
equal  to  zero.  This  approach  results  in  a  high  accuracy  interpolation  between  points  in 
wavenumber  using  the  discrete  Fourier  transform.  The  data  are  transformed  into  the  wavenumber 
domain  using  the  equations 
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=  =  Y,w{XnMxn,a)o)e 


2jt 


(3.8a) 


n=0 


and 


= 


N  \27t 


2% 


(3.8b) 


M  +  N)  L  {M-¥N)K^ 
where  w(x„)  is  the  spatial  weighting  function  (window)  of  the  discrete  Fourier  transform  evaluated 
at  position  and  Q)q  represents  a  fixed  frequency.  Changing  the  window  function  allows  data  to 
be  processed  with  different  sidelobe  levels  and  different  mainlobe  widths.  The  weighting  function 
was  normalized  such  that 


N+M-l 
]^w(x„)  =  l  . 
n=0 


(3.9) 


3.5.  FREQUENCY-WAVENUMBER  DOMAIN  MEASUREMENTS 

The  discrete  Fourier  transform  described  previously  was  applied  to  the  spatial  domain 
measurements.  Two  aperture  sizes  were  processed:  one  with  n  =  15  channels  and  a  separation 
distance  of  0.00635  m  (0.25  inch)  and  another  with  8  channels  and  a  separation  distance  of 
0.0127  m  (0.50  inch).  Figures  3.8  to  3.27  illustrate  the  wavenumber  cuts  of  the  transformed 
spatial  domain  data  starting  at  100  Hz  and  ending  at  480  Hz  in  20-Hz  increments.  These  data  are 
displayed  as  a  function  of  wavenumber  using  a  uniform  (rectangular)  window  with  each  individual 
weighting  function  equal  to  0.005  (=  1/{N+M))  and  a  zero  pad  of  185  (zeros)  to  15  (actual  data 
points)  for  the  15-channel  aperture  and  a  zero  pad  of  192  to  8  for  the  8-channel  aperture  (N+M  = 
200).  The  wave  energy  and  wave  speeds  are  identified  on  the  plots.  The  wave  speeds  were 
determined  using  the  equation 


c 


(3.10) 
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where  c  is  the  apparent  wave  speed  (m/s)  sensed  by  the  array,  a  term  also  frequently  referred  to  as 
the  phase  velocity.  Figure  3.28  is  a  plot  of  the  theoretical  shear  wave  speed  of  Hexcel  195-RE 
urethane  and  the  measured  apparent  wave  speeds  of  the  array  versus  frequency.  The  theoretical 
shear  wave  speed  was  determined  using  (Verburg,  1983) 

(3.11) 

where  is  the  shear  wave  speed,  G  is  the  shear  modulus,  and  p  is  the  density.  The  area  where 
the  theoretical  wave  speed  is  decreasing  with  respect  to  frequency  is  shown  (dashed  line)  because 
this  is  inconsistent  behavior  for  the  material  according  to  many  material  laws.  Note,  however,  that 
the  measured  shear  wave  speed  is  decreasing  (with  respect  to  frequency)  near  the  frequencies 
where  the  material  measurements  predict  a  decreasing  shear  wave  speed.  Figure  3.29  shows  the 
wavenumber-frequency  surface  of  the  measurements. 

From  figure  3.1,  it  is  observed  that  turbulent  flow  induces  a  maximum  spectrum  level  in 
the  100-  to  125-Hz  region.  An  examination  of  figures  3.2,  3.3,  and  3.4  indicates  that  there  is  a 
relative  minimum  in  the  amplitude  transfer  function  between  the  reference  sensor  and  the  locations 
at  1.09,  2.75,  and  4.50  inches.  A  likely  cause  of  these  local  minima  in  the  transfer  function  is 
coherent  destructive  interference  between  different  propagation  modes  (paths). 

Consider  the  case  of  two  coherent  modes  contributing  to  the  waveform, /„(0,  observed  at 

the  nth  sensor  location.  This  waveform  can  be  expressed  as 

~  —  T^l)  ~  ^n2)  ^n(^)  ’  (3.12) 

where /(r)  is  the  source  waveform  at  sensor  1,  t„j  is  the  effective  propagation  time  relative  to 
sensor  1  from  the  source  to  the  nth  sensor  supported  by  the  first  propagation  mode  (path),  and  ^n2 

is  the  effective  propagation  time  relative  to  sensor  1  from  the  source  to  the  nth  sensor  supported  by 
the  second  mode  (path).  The  parameter  a^p  is  the  relative  attenuation  of  amplitude  for  the  p\h 
mode  relative  to  the  first  mode  referenced  to  the  nth  sensor.  The  waveform  v^(r)  is  the  additive 
noise  at  the  nth  sensor.  Equation  (3.12)  has  the  Fourier  transform 
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Xnico)  = 


CCfii  +  (Xfi2^ 


(3.13) 


where  F{oS)  is  the  Fourier  transform  of  fit)  and  A„  12  =  '^nl  ~  '^nl-  APSD  for  the  nth  sensor  is 

(3.14a) 


J?„(0)  =  |X„(< 


or 


Rjinico)  =  \F{co)f[ali  +  al2  +  2a„ia„2  cos©A„,i2) +|V„(fi))|' 
The  transfer  function  from  sensor  1  to  sensor  n  is 


Hinico) 


Xfi(0)Xn{(0) 


Ril(co) 


which  is  also 


Hi„(ty)=- 


OCfii  +  CCji2^ 


J(oAn,l2 


-mni-Tnl) 


OC11  +  ai2e^^^^’^^]  + 

.  .0 

/y  ^ry  ^ 

J 

J 

(3.14b) 


(3.15a) 


(3.15b) 


where  \F{(o)f' /^\{(0)f  is  the  SNR  at  frequency  (O  (rad/s).  At  high  SNR,  equation  (3.15b)  gives 
the  squared  amplitude  of  the  transfer  function  as 

and  shows  the  linear  dependence  of  phase  with  respect  to  frequency.  It  is  argued  that  because 
sensor  1  is  located  on  the  urethane  block  surface  at  the  point  closest  to  the  blockage,  there  is 
predominantly  only  a  single  path  from  the  source  to  sensor  1,  as  shown  in  figure  3.30.  Hence, 
Til  =  '^12  Ai  12  =  Til  -  112  =  0.  Thi^  argument  is  supported,  at  least  above  60  Hz,  by  the 

fact  that  there  are  no  nulls  in  the  APSD  of  figure  3.1  that  would  indicate  intermodal  destructive 
interference.  For  sensor  locations  n  =  2,  3,  4  .  . . ,  equation  (3.16)  predicts  a  local  minimum  at 
©jninAn  12  =  (2^  + 1)^-  For  a  mean  propagation  velocity  of  c  =  13  m/s,  an  effective  propagation 

path  length  difference,  12  ,  between  modes  1  and  2  is  seen  to  be 


«nl  +  cell  +  2«niQ:n2  COSt»A„j2 
ail  +  ®12  +  2aiiai2  costoAi  12 


(3.16) 
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(3-17) 

From  figure  3.3,  with  a  sensor  separation  of  2.75  inches  (6.99  cm),  the  local  minimum  in 
the  transfer  function  at  approximately  312  Hz  predicts  a  k  =  0  mode  interference  effective  path 
length  difference  of 

5„42(0)  =  0-82  inch  (2.1  cm)  .  (3.18) 

Similarly,  for  the  weaker  SNR  null  at  360  Hz  in  figure  3.4  (separation  4.5  inches  (1 1.43  cm)),  an 
effective  path  length  difference  of  0.75  inch  (1.9  cm)  is  predicted.  For  the  assumed  two  mode 
model,  the  effective  path  length  differences  calculated  above  for  the  two  propagation  modes  are 
entirely  consistent  with  the  test  geometry.  That  is,  the  effective  path  length  differences  are  less 
than  any  direct  path  propagation  distance.  In  fact,  the  calculated  effective  path  differences  are 
supportive  of  a  model  that  includes  a  direct  path  and  a  second  path,  perhaps  a  surface  wave,  from 
the  source  to  the  sensors.  However,  a  more  exact  understanding  of  the  actual  spatial  extent  of  the 
turbulence-induced  excitation  source  and  the  full  solution  of  the  frequency  domain  (Helmholtz) 
wave  equation  is  clearly  needed  to  substantiate  the  very  basic  approach  taken  above. 
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Figure  3.1 1.  Wavenumber  Cut  at  160  Hertz 
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Figure  3.29.  Wavenumber-Frequency  Surface  of  Measured  Data 
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Figure  3.30.  Geometry  of  Four  Sensor  Locations  Used  to  Measure  Transfer 
Functions  Between  Sensor  Locations  (1  to  2),  (1  to  3),  and  (1  to  4).  Energy 
Propagation  From  the  Source  to  Sensor  Location  1  is  Assumed  To  Be  Single 

Mode,  Direct  Path  (pi). 


3.6.  SUMMARY 

The  propagating  energy  produced  by  the  blockages  in  the  tubing  can  be  measured  with 
sensors  applied  to  the  surface  of  the  urethane  solid.  Transferring  the  data  into  the  wavenumber 
domain  provides  a  method  to  identify  wave  speeds  of  coherent  energy  in  the  solid  block. 
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CHAPTER  4.  PASSIVE  IMAGING  OF  A  SIMULATED  EXTENDED 

NEARFIELD  RADIATOR 

4.1.  INTRODUCTION 

The  array-source  geometry  for  the  simulation  of  an  extended  nearfield  source  is  depicted  in 
figure  4.1.  The  eight-sensor  line  array  is  located  on  the  x-axis,  and  all  dimensions  are  in 
centimeters.  A  radiating  source  is  simulated  using  a  total  of  64  point  sources  distributed  in  the 
cylinder.  The  sources  were  equally  spaced  in  each  direction:  radius,  length  of  the  cylinder,  and 
circular  cross  section.  Each  point  source  signal  and  received  sensor  noise  were  simulated  using 
independent  normal  Gaussian  random  numbers  and  a  sampling  rate  of  2000  Hz.  The  received 
signal  at  each  passive  sensor  is  the  sum  of  signals  from  each  of  the  point  sources  and  noise. 
Linear  time  delay  interpolation  and  a  wave  speed  of  10  m/s  were  used  to  calculate  the  delay  from 
each  source  to  each  sensor.  From  the  geometry,  it  is  clear  that  the  effective  source  is  extended  and 
nearfield.  The  received  sensor  signal  was  filtered  using  a  Butterworth  filter  of  order  14,  with  a 
bandpass  from  100  to  350  Hz  and  a  stopband  from  50  to  400  Hz. 

4.2.  MUSIC  AND  MAGDR  SIGNAL  PROCESSING 

Eigenvectors  and  eigenvalues  of  the  sensor-to-sensor  cross  spectral  density  matrix  (CSDM) 
were  obtained  using  the  MATLAB  singular  value  decomposition  (SVD)  method.  The  eight-by¬ 
eight  CSDM  was  estimated  at  250  Hz  by  using  a  256-point  block  fast  Fourier  transform  (FFT)  and 
a  40-snapshot  (5.12  sec)  average.  A  Hamming  window  was  applied  to  the  input  sample  time  data 
before  the  FFT  transform.  The  eigenvalue  amplitudes  for  sensor  output  SNR  ratios  of  12  and  3 
dB  are  plotted  in  figures  4.2  and  4.3,  respectively.  The  maximum  array  gain  distortionless 
(MAGDR)  beamformer  and  MUSIC  direction  finder  (DF)  responses  were  calculated  for  SNR 
levels  of  12  and  3  dB.  For  the  MUSIC  response,  the  eigenvectors  corresponding  to  the  three 
dominant  eigenvalues  were  used  as  the  source  vector  subspace.  The  MAGDR  beamformer  and 
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MUSIC  DF  normalized  response  images,  viewed  at  a  vertical  angle  of  30  degrees,  are  presented  in 
figures  4.4-4.15  with  a  spatial  resolution  of  0.25  cm  in  each  dimension. 

Figures  4.4,  4.8,  and  4.12  are  the  MUSIC  DF  response  functions  in  the  x-y,  y-z,  and  x-z 
planes,  respectively,  for  a  sensor  output  SNR  of  12  dB.  Figures  4.5,  4.9,  and  4.13  are  the 
MAGDR  beamformer  response  functions  in  the  x-y,  y-z,  and  x-z  planes,  respectively,  for  a  sensor 
output  SNR  of  12  dB.  Figures  4.6,  4.10,  and  4.14  are  the  MUSIC  DF  response  functions  in  the 
x-y,  y-z,  and  x-z  planes,  respectively,  for  a  sensor  output  SNR  of  3  dB.  Figures  4.7,  4.11,  and 
4.15  are  the  MAGDR  beamformer  response  functions  in  the  x-y,  y-z,  and  x-z  planes,  respectively, 
for  a  sensor  output  SNR  of  3  dB.  All  the  response  planar  cuts  were  chosen  to  intersect  the  source 
volume. 

The  following  conclusions  can  be  drawn  from  the  response  plots: 

•  The  largest  eigenvalue  becomes  less  dominant  as  the  sensor  output  SNR  decreases. 

•  A  cusp  of  ambiguity  for  response  plots  in  the  x-z  plane  occurs  because  of  the  conical 
response  pattern  for  the  linear  sensor  array. 

•  The  y-z  response  plot  has  some  feature  in  the  z-direction  due  to  the  nearfield.  This 
perspective  is  the  most  informative  because  it  gives  the  best  indicator  of  distance  (2  cm)  to  the 
source  for  z  =  0  and  allows  an  approximation  of  10  log  8  =  9  dB  for  the  SNR  improvement  due  to 
beamforming  (see  figure  4.12). 

•  For  the  3-dB  SNR  case,  the  spatial  response  results  in  high  sidelobes  and  ambiguity 
(compare  figure  4. 14  to  figure  4.12). 
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Figure  4.4.  Normalized  MUSIC  Response  in  the  x-y  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  0.25,  y  =  1.25)  and  0.103  (x  =  1.25,  y  =  10.0) 
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Figure  4.5.  Normalized  MAGDR  Response  in  the  x-y  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  0.25,  y  =  1.5)  and  0.3536  (x  =  5.0,  y  =  2.25) 
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Figure  4.6.  Normalized  MUSIC  Response  in  the  x-y  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  4.5,  y  =  0.0)  and  0.2395  (x  =  3.5,  y  =  0.0) 
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Figure  4.7.  Normalized  MAGDR  Response  in  the  x-y  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  2.25,  y  =  1.0)  and  0.477  (x  =  5.0,  y  =  9.0) 
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Figure  4.8.  Normalized  MUSIC  Response  in  the  x-z  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  -0.25,  z  =  0.0)  and  0.12  (x  =  5.0,  z  =  -3.75) 
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Figure  4.9.  Normalized  MAGDR  Response  in  the  x-z  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  -0.25,  z  =  -0.5)  and  0.362  (x  =  5.0,  z  =  -1.0) 
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Figure  4.10.  Normalized  MUSIC  Response  in  the  x-z  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  2.25,  z  =  0.0)  and  0.431  (x  =  -1.0,  z  =  -5.0) 
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Figure  4. 1 1 .  Normalized  M  AGDR  Response  in  the  x-z  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (x  =  2.25,  z  =  0.0)  and  0.5203  (x  =  5.0,  z  =  -5.0) 
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Figure  4. 12.  Normalized  MUSIC  Response  in  the  y-z  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (y  =  0.5,  z  =  -1.5)  and  0.113  (y  =  10.0,  z  =  -5.0) 
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Figure  4.13.  Normalized  MAGDR  Response  in  the  y-z  Plane  for  a  Sensor  Output 
SNR  of  12  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (y  =  0.25,  z  =  -1.7)  and  0.4557  (y  =  6.25,  z  =  -4.5) 
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Figure  4.14.  Normalized  MUSIC  Response  in  the  y-z  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (y  =  0.5,  z  =  -0.75)  and  0.6665  (y  =  1.25,  z  =  -4.0) 
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Figure  4.15.  Normalized  MAGDR  Response  in  the  y-z  Plane  for  a  Sensor  Output 
SNR  of  3  dB  at  250  Hz.  Maximum  and  Minimum  Response  Values 
Are  1.0  (y  =  0.0,  z  =  -1.25)  and  0.8292  (y  =  0.75,  z  =  -4.5) 
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CHAPTER  5.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  results  presented  in  this  report  show  that  the  turbulent  energy  induced  into  a  solid 
viscoelastic  medium  can  be  detected  and  that  the  energy  can  be  resolved  using  a  conventional 
beamformer.  Additionally,  simulation-based  model  predictions  show  that  this  type  of  source  has 
significant  potential  for  being  localized  and  imaged.  The  laboratory  experimental  work  should  be 
expanded  to  measurements  taken  with  background  noises,  nonhomogeneous  media,  advanced 
beamformers,  and  a  sensor  array  rather  than  a  series  of  two-point  measurements.  This  approach 
will  provide  a  better  understanding  of  such  issues  as  SNR,  detection  threshold,  and  human 
implementation. 

In  order  for  this  type  of  detection,  localization,  and  imaging  to  be  executed  in  humans,  the 
following  technical  issues  wiU  have  to  be  resolved: 

(a)  Transducer  Phase  Accuracy  Requirements  -  Array  transducer  location  accuracy  is 
absolutely  essential.  Such  array  geometry  accuracy  will  likely  dictate  either  a  rigid  or  semirigid 
frame  to  which  the  transducers  can  be  attached.  It  is  envisioned  that  this  frame  can  be  suspended 
in  a  completely  reusable  exterior  air  bladder,  which  could  be  inflated,  thereby  forcing  the  sensor 
array  to  conform  tightly  to  the  patient’s  skin  in  order  to  maximize  energy  coupling  and  reduce 
impedance  mismatch.  Acoustic  absorption  at  the  rear  of  the  bladder  must  be  maximized  to 
minimize  reflected  energy  interference.  This  issue  must  be  given  the  primary  attention. 

(b)  Signal  and  Masking  Source  Acoustic  Power  Levels  Relative  to  the  Statistically 
Uncorrelated  Component  of  the  Transducer  Output  -  The  relative  levels  of  the  various  signal 
sources  and  system  noise  with  respect  to  the  vibration  sensing  transducer  pickup  noise  component, 
which  is  statistically  uncorrelated  between  transducers,  are  major  determinants  of  the  extent  to 
which  desired  waveforms  can  be  extracted  from  the  background  noise  and  subsequently  resolved 
spatially  from  other  concurrent  interfering  sources. 

(c)  Affected  Artery  Segments  and  Sources  of  Acoustic  Masking  Energy  -  The  limits  of 
spatial  resolution  and  the  required  complexity  of  the  digital  signal  processing  (DSP)  hardware  are 
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determined  by  the  angular  extent  of  the  particular  distribution  of  the  sources  of  energy.  In 
addition,  the  spatial  coherence  of  the  extended  source  is  significant.  For  example,  the  nonlinear 
(chaotic)  nature  of  turbulent  flow  either  through  a  stenotic  valve  or  regurgitant  valve  would  tend  to 
be  spatially  incoherent,  whereas  the  pressure  field  generated  by  the  synchronized  movement  of 
points  on  a  vessel  wall,  which  is  being  excited  to  resonance  by  turbulent  internal  flow,  may  be 
spatially  coherent. 

(d)  Temporal  Extent  for  Which  Cycle-to-Cycle  Acoustic  Energy  Can  Be  Coherently 
Integrated  -  Energy  field  temporal  coherence  can  be  used  to  enhance  high  performance  modem 
DSP  algorithm  capability  if  the  received  space-time  field  wave  form  can  be  synchronized  from  one 
heart  cycle  to  the  next. 

(e)  Cost-Effective  Hardware-Software  Realization  -  Modem  high  performance  DSP 
algorithms  are  notably  very  computationally  intense  with  regard  to  matrix  (linear)  algebra 
requirements.  While  arithmetically  intensive,  such  matrix- vector  computations  are  very  amenable 
to  efficient  implementation  on  parallel  architecture  computing  devices. 

(f)  Homogeneity  of  the  Sound  Propagation  Medium  -  The  capability  to  realize  both  a  high 
SNR  gain  and  high  resolution  is  limited  by  the  amount  of  uncertainty  in  the  propagation  time 
assumed  from  a  point  of  focus  either  on  or  in  the  heart  to  each  of  the  auscultation  points.  This 
propagation  time  is  dependent  on  the  degree  of  uniformity  in  material  density  along  the  propagation 
path.  For  example,  fatty  tissue  would  presumably  be  less  dense  than  fluid  and  would  therefore 
exhibit  a  slower  propagation  velocity.  Moreover,  the  absorption  and  scattering  of  energy  may 
complicate  the  vibration  field  beyond  what  can  be  easily  modeled  as  straight-line  propagation.  This 
problem  may  require  a  complete  and  accurate  solution  to  the  wave  equation  within  the  torso,  using 
the  appropriate  boundary  conditions. 

It  should  be  explicitly  noted  that  the  energy  sources  are  located  in  the  very  nearfield  of  the 
auscultation  transducer  array.  Moreover,  because  of  the  proximity  of  source  locations  and  the 
transducer  array,  the  traditional  radar/sonar  assumption  of  a  spatially  stationary  signal  field  with 
uniform  spatial  placement  of  the  transducers  in  the  array  may  be  violated.  This  aspect  of  the 
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problem  eliminates  from  consideration  the  use  of  high  performance  locator  algorithms,  such  as 
autoregressive  (AR),  autoregressive  moving  average  (ARMA)  (Akay,  1990)  and  ESPRIT 
(Owsley,  1985;  Owsley,  1992).  Specifically,  the  reduced  rank  minimum  variance  distortionless 
response  (MVDR)  beamfoming  algorithm  is  a  prime  candidate  for  waveform  extraction  and  the 
enhanced  MVDR  algorithm  is  a  valid  choice  for  the  location  of  lesion-induced  resonant  points  in 
the  coronary  artery  mesh  (Owsley,  1985;  Owsley,  1992). 
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